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We show that single-particle energies in doubly magic nuclei depend almost linearly on the cou- 
pling constants of the nuclear energy density functional. Therefore, they can be very well char- 
acterized by the linear regression coefficients, which we calculate for the coupling constants of the 
standard Skyrme functional. We then use these regression coefficients to refit the coupling constants 
to experimental values of single-particle energies. We show that the obtained rms deviations from 
experimental data are still quite large, of the order of 1.1 MeV. This suggests that the current stan- 
dard form of the Skyrme functional cannot ensure spectroscopic-quality description of single-particle 
energies, and that extensions of this form are very much required. 
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I. INTRODUCTION 



The energy-density-functional (EDF) methods are re- 
cently intensely studied in nuclear physics to gain more 
precise and predictive description of stable and exotic 
nuclei within a unified framework. Such methods can be 
considered as realizations of the nuclear effective theory 
at low energies [l[ and are motivated by successful ap- 
plications of the density functional theory in electronic 
systems, see, e.g., Ref. @- In nuclear physics, they 
have been studied since early 1970s under the names of 
Hartree-Fock or Hartree-Fock-Bogolyubov methods [|[, 
but in fact they rather correspond to the Kohn-Sham 
[H type of approaches, aiming at describing correlated 
fermion systems by their one-body (local or non-local) 
densities only. 

In the present study we focus on the Skyrme-type 
EDFs, by which one assumes that the ground-state en- 
ergy of a given nucleus is given by an integral of a local 
energy density. As discussed in Ref. ||, our main goal is 
to look for the spectroscopic-quality EDF, which would 
correctly describe positions of single-particle (s.p.) en- 
ergies across the nuclear mass table. To this end, here 
we analyze dependence of the s.p. energies on the EDF 
coupling constants and attempt answering the question 
on whether the current parametrizations of the EDF are 
rich enough to describe experimental data with reason- 
able precision. In the present work we look only at bare 
s.p. energies in doubly magic nuclei, because the polar- 
ization effects, which may affect detailed values of s.p. 
energies, are estimated to be significantly smaller than 
discrepancies with experimental data [5j. 

The paper is organized in the following way. In Sec. |TT] 
we briefly introduce necessary definitions and present our 
method of analysis. Three subsections of Sec. IIIII present 
our results on s.p. energies in function of coupling con- 
stants, regression coefficients, and fits to experimental 
data. Conclusions of our work are presented in Sec. IIVI 



II. METHOD 

In the present study we consider EDF in the form given 
in Refs. [1,0], 



£ = 



J d 3 rH(r), 



(1) 



where the energy density Tt(r) can be represented as a 
sum of the kinetic energy and of the potential-energy 
isoscalar (t = 0) and isovector (t = 1) terms, 



H(r) = —T + H (r)+H 1 (r), 
Zm 



(2) 



which for the time-reversal and spherical symmetries im- 
posed read: 

H t {r)=C p tP 2 t +C^p t &p t +C T t p t T t + \C J t J 2 t +C7 J p t V-J t . 

Standard definitions of the local densities pt, r t , and Jt 
were given in Refs. [7J, |8( and are not repeated here. Fol- 
lowing the parametrization used for the Skyrme forces, 
we assume the dependence of the coupling parameters 
C t p on the isoscalar density po as: 

qf = + (4) 

Similarly as in Ref. |H, we note here that EDF of Eq. (J) 
depends linearly on twelve coupling constants, 



r<f> 
°to> 
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(5) 



for t = and 1. Therefore, due to the Hellmann-Feynman 
theorem derivatives of the total energy E with re- 
spect to the coupling constants are given by space inte- 
grals of densities appearing in Eq. © 0. 

Variation of EDF in Eq. (JXJ) with respect to the Kohn- 
Sham orbitals 4>i(r), which define the local densities, 
gives the standard eigenequation, 



h<t>i = e l cf) l 



(6) 



2 



where the Kohn-Sham one-body Hamiltonian h is ob- 
tained as a functional derivative of EDF with respect 
to densities, see, e.g., Refs.@, H, EH- ^ i s obvious that 
h also depends linearly on the coupling constants, and 
therefore, again due to the Hcllmann-Fcynman theorem, 
derivatives of the s.p. energies £j with respect to the cou- 
pling constants are also given by space integrals of densi- 
ties. As a consequence, we can expect that these deriva- 
tives are generic quantities weakly dependent on a par- 
ticular parametrization of the functional, provided the 
functional has been adjusted to data and the underlying 
densities are correct. 

The aim of the present work is not only to determine 
these derivatives but also to estimate to which extent 
they are generic. Of course, for a given set of coupling 
constants, they can be calculated from the Hellmann- 
Feynman theorem. Equivalently, their determination 
may rely on numerically calculating functions ej(C m ), 
where index m = 1, ... 12 is assumed to enumerate the 
twelve coupling constants ([5]). One simple option would 
be to calculate them from the finite-difference formula, 
for example as 



(Cm) 



j(Cm) — e i(C m ) 



(7) 



Cm Cm 

with values of C+ and C~ suitably close to C^. Had 
the functions ei(C m ) been exactly linear, the Hellmann- 
Feynman and finite-difference methods would have given 
exactly the same answers. If we aim at determining the 
degree to which they are not linear we have to proceed 
in another way. 

To this end, the method of choice is the linear regres- 
sion analysis [l2j , which makes a hypothesis of linearity 
and tests it by determining the regression coefficients and 
their standard deviations. In our case, we write the ex- 
pression 

e i {C m +d k m )=I im d k m + I° i +r k i , (8) 

where Ci(C^ n + d^) are the s.p. energies calculated self- 
consistently for coupling constants + d^ that differ 
from by suitably small numbers d^. Index k enu- 
merates different choices of these small numbers (which 
sample the vicinity of C^J and r\ are the residuals be- 
tween the self-consistent results and linear approxima- 
tion given by the regression coefficients Ii m and e°. The 
regression method minimizes the residuals by adjusting 
the regression coefficients and determines their standard 
deviations. It is obvious that regression coefficients Ii m 
constitute estimates of derivatives ([7]), while if are very 
close to €i(Cj n ). 

The results below were obtained by using d^ = ±dC^, 
where d = 0.001 to 0.005, depending on the Skyrme 
force and nucleus, i.e., each of the twelve coupling con- 
stants was raised or lowered by the same percent frac- 
tion. (In some cases, for vanishing coupling constants 
C^j, appropriate absolute shifts were used.) As a result, 
in each case the regression analysis was done by using 
k = 1,.. . ,2 12 =4096 samples. 
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FIG. 1: (Color online) Neutron (left) and proton (right) s.p. 
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levels in Ca in function of 5C- 

C£ (SLy5)=-931.294MeVfm 3 (as well as other, unchanged, 
coupling constants) correspond to the Skyrme force SLy5. 
Symbols listed in the legend from left to right correspond to 
increasing s.p. energies. The range of the d parameter used 
in the regression analysis is also indicated. 



III. RESULTS 

In the present work, the effect of the EDF coupling 
constants on s.p. energies was studied in doubly magic 
nuclei. Calculations were performed by using the code 
HFBRAD 13], which solves the self-consistent equa- 
tions with spherical symmetry assumed. Nuclei consid- 
ered in our calculations were 4 He, 16 , 40 - 48 Ca, 48 ' 56 Ni, 
100 - 132 Sn, and 208 Pb. The coupling constants C^ © cor- 
respondin g to the Skyrme functionals SLy5 [l4| , SkP [l5| , 
and SkO' [id ] were used. In order to consistently calcu- 
late the energies of the occupied and unoccupied states, 
the program was run in the HFB mode. Even if the pair- 
ing correlations vanish in doubly-magic nuclei, a residual 
non-zero pairing was kept and equivalent s.p. energies 
[1511 were determined. 



A. Single-particle energies in function of coupling 
constants 

In Figs. [THS] we show the s.p. energies in 40 Ca as func- 
tions of five isoscalar coupling constants ([5]). Results 
were obtained by varying one coupling constant of the 
SLy5 functional and by keeping all the remaining ones 
at their SLy5 values. Ranges of variation of the coupling 
constants were chosen in a maximum possible way, i.e., 
up to the values where changes of level ordering or lev- 
els becoming unbound precluded obtaining meaningful 
solutions. Of course, physical values of these coupling- 
constants are fairly well fixed by adjustments to empiri- 
cal data, so, in practice, physically acceptable variations 
must not be so large. Therefore, values of parameters 
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FIG. 2: (Color online) Same as in Fig. Q] but for 5C=C { 
C*£ D (SLy5), where C^ D (SLy5)=859.813 MeV fm 3(1+a) 
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FIG. 3: (Color online) Same as in Fig. [T] but for 5C=C5 
CJ(SLy5), where C5 (SLy5)=56 MeV fm 5 . 
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FIG. 4: (Color online) Same as in Fig. Q] but for SC=Co" 
C ( f p (SLy5), where C Ap (SLy5)=-76.793 MeV fm 5 . 



FIG. 5: (Color online) Same as in Fig. Q] but for 5C=C% 3 
C VJ (SLy5), where C VJ (SLy5) = -93.75 MeV fm 5 . 



d used in the regression analysis below, also shown in 
Figs. HHSl were quite small. Positive s.p. energies are 
shown only to indicate in which regions of parameters the 
levels become unbound - precise values of these s.p. en- 
ergies are only loosely related to positions of resonances. 

In Figs. HH5] one can clearly see that there is only a 
rather small overall non-linearity of the s.p. energies as 
functions of the coupling constants. With decreasing val- 
ues of coupling constants Cq and Cp D , the s.p. potentials 
become deeper and thus the values of s.p. energies uni- 
formly decrease (Figs. []] and [5]). There are only very 
small differences in the dependencies induced by vary- 
ing coupling constants Cq and Cq D . Somewhat larger 
changes in relative positions of s.p. levels are induced 
by varying coupling constant Cq (Fig. [3]). Decreasing 
values of this coupling constant induce deeper s.p. po- 
tentials and larger values of the effective mass. Coupling 
constant C p clearly influences the surface properties of 
the s.p. potentials (Fig. [4j) by changing relative positions 
of the low-£ and high-£ levels. Finally, the SO coupling 
constant Cq' 3 very linearly changes the SO splitting of 
levels shown in Fig. Dependence of levels on the ten- 
sor coupling constant Cg is very weak in a spin-saturated 
nucleus 40 Ca, and, therefore, it is not shown. 



B. Regression coefficients 

It is clear that results shown in Figs. [T]-[5] can only 
tentatively indicate the type of dependence of the s.p. 
energies on coupling constants. In order to have a better 
handle on this dependence in a multi-dimensional space 
of parameters, it is more efficient to directly consider the 
regression coefficients discussed in Sec. |TTJ For a better 
separation of central and SO effects, the regression anal- 
ysis was performed independently for centroids and SO 
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FIG. 6: (Color online) Bar charts of regression coefficients 
Iim, Eq. (|8]), for neutron (left) and proton (right) centroids 
(top) and SO splittings (bottom), determined for the 2d or- 
bital and SLy5 Skyrme functional. The order of bars is the 
same as the order of nuclei in the legend on top of the Figure. 
In order to better compare results for the centroids and SO 
splittings, regression coefficients are shown for inverted (neg- 
ative) SO splittings. Units of regression coefficients 7j m are 
equal to a MeV divided by the units of coupling constants, 
which are given in captions of Figs. [JJ-[5] 



splittings of the s.p. levels, 
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(9) 
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respectively. Below we present selected results illustrat- 
ing the main conclusions of the present study. 

Regression coefficients obtained for the centroids and 
SO splittings of the Id levels in five light nuclei are shown 
in Fig. El As expected, in the N — Z nuclei 16 O, 40 Ca, 
and 56 Ni, the s.p. energies depend only on the isoscalar 
coupling constants. Only the centroids in N ^ Z nu- 
clei 48 Ca and 48 Ni weakly depend on the isovector cou- 
pling constants Cf , Cf D , and CJ", and the SO splittings 
in these nuclei weakly depend on the isovector tensor 
coupling constant Cf. In all cases, dependence on the 
isovector SO coupling constant C^ J is extremely weak. 
This illustrates the fact that the isospin excess in 48 Ca 
and 48 Ni is still not large enough for a pronounced de- 
pendence on the isovector coupling constants, and thus 
these coupling constants can be reasonably fixed only by 
going to even more exotic nuclei. 

Centroids of the Id levels depend predominantly on 
four isoscalar coupling constants, Cq , Cq D , CJ, and Cq P . 
These coupling constants act on positions of centroids in 
a fairly similar way, i.e., by looking at positions of Id 
orbitals in different nuclei one is not able to distinguish 
between central, effective mass, and surface effects, or 
differentiate between the density-dependent and density- 
independent coupling constants. 

In all studied nuclei we obtain a standard uniform 
dependence of the SO splittings on the isoscalar SO 



FIG. 7: (Color online) Same as in Fig.[6]but for the regression 
coefficients l' im , Eq. (llip . 



coupling constant . A somewhat weaker, but non- 
negligible, dependence on the isoscalar tensor coupling 
constant Cq is also clearly seen in spin-non-saturated nu- 
clei 48 Ca, 48 Ni, and 56 Ni. An unexpected result of our 
study is the fact that the SO splittings also very strongly 
depend on the central, Cq and Cq D , effective-mass, Cq , 
and surface, C^ p coupling constants. These dependen- 
cies were not explicitly recognized in the literature, be- 
cause most often the SO splittings were adjusted at the 
end of the fitting protocol, i.e., for fixed values of other 
coupling constants. It is clear, however, that a global ad- 
justment will create specific inter-relations between the 
SO and tensor coupling constants on one side and all the 
other ones on the other side. Relations of these kind 
between the SO and effective-mass effects were recently 
analyzed in Ref. p7| . 

In order to elucidate origins of the obtained depen- 
dencies of the SO splittings on the coupling constants, 
in Fig. [7] we show results of the regression analysis per- 
formed for the s.p. energies relative to those of the Is 
orbital, i.e., in analogy with Eq. ©, 



i(Cm + dm) 
■(^m + ^m) 



r d k 

irn 77 



(ll) 



In this way, we remove the overall effects of scaling of the 
s.p. spectra related to changing the depth of the central 
potential and/or effective mass. 

One can see that, indeed, in this changed representa- 
tion, dependence of relative centroids on coupling con- 
stants becomes significantly weaker. Also dependence of 
the relative SO splittings on the central, Cg and Cq D , 
and surface, Cq P , coupling constants becomes less pro- 
nounced. Nevertheless, the relative SO splittings still 
quite strongly depend on the effective-mass coupling con- 
stant, Cq, i.e., this dependence cannot be attributed to 
a simple scaling of the s.p. spectrum. 

This fact is already clearly visible in Fig. 02 where 
the s.p. energies of favored SO partners, lds/2 and lfj/2. 
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FIG. 8: (Color online) Derivatives of the particle (in fm -3 ), 
kinetic (in fm -5 ), and spin-orbit (in fm -4 ) neutron densities 
with respect to the coupling constant Cq (in MeVfm 5 ). The 
inset shows the densities themselves. Calculations were per- 
formed in 40 Ca for the SLy5 functional. 
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FIG. 9: (Color online) Same as in Fig. [8] but for the central 
(in MeV), spin-orbit (in MeVfm), and kinetic (in MeVfm 2 ) 
potentials. 



clearly follow that of the 2si 
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unfavored SO partners, ld 3 / 2 and lf5/2, 



level, while those of the 
behave quite dif- 
ferently. 

The influence of the isoscalar effective-mass coupling 
constant on the SO splittings can be understood by in- 
specting the self-consistent densities and potentials. In 
Fig. [8]we show derivatives of the kinetic (a), particle (b), 
and SO (c) neutron densities in 40 Ca with respect to the 
coupling constant CJ. For completeness, the inset shows 
the total densities. The derivatives were determined from 
differences of densities calculated self-consistently at two 
values close to the SLy5 value of Cq. One can see that 
an increase in Cq not only lowers the kinetic density r n , 
but also significantly lowers the particle density p n in the 
center of the nucleus and moves the particles to the sur- 
face. As a result, it induces a smaller gradient of the 
particle density at the surface, and hence a weaker SO 
potential, as shown in Fig. [9l and as a consequence - a 
smaller SO splitting. 

One can say that the inter-relations between the 
effective-mass and SO effects stem from the fact that 
EDF of Eq. ([3]) contains the effective-mass term, CqPqTq, 
that depends on the product of the particle and kinetic 
densities. Therefore, changes in the effective-mass cou- 
pling constant induce changes in both these densities, 
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FIG. 10: (Color online) Same as in Fig.[6]but for the If orbital. 
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FIG. 11: (Color online) Same as in Fig. [6] but for the lg 
orbital. 



and, in particular, they induce changes in the gradient 
of the particle density that affect the SO splittings. The 
role of the effective-mass term is clearly not restricted to 
changing the overall density of levels, as is the case in a 
trivial case of the infinite square-well potential, cf. also 
discussion in Ref. [T^. 

Conclusions drawn from the results obtained for the 
Id levels are further corroborated by those for the If, 
lg, and lh levels, which are shown in Figs. [Hjl [Til an d 
Q2J respectively. One can see that a tangible influence of 
the isovector coupling constants on centroids is obtained 
only in 132 Sn and 208 Pb, i.e., in systems with a signifi- 
cant neutron excess. But even there the influence of the 
isovector coupling constants on the SO splittings is quite 
weak. Again we see that adjustment of the isovector cou- 
pling constants to empirical data may require studying 
systems much further away from stability. 

A universality of the regression coefficients Ii m , 
Eq. ([H]), is illustrated in Figs. [T3] and [TJJ in which we 
present results obtained for the SkP and SkO' Skyrme 
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FIG. 12: (Color online) Same as in Fig. [6] but for the lh 
orbital. 



0.3 



0.2 



0.0 



16„ 40„ 48„ 48... 56... 

i O ■ Ca □ Ca ■ Ni ■ Ni 



5b n 1 
Pi 

0.0 
-0.1 



neutrons (l/2)(ld 5/2 +ld 3/2 ) 



protons (l/2)(ld g/2 +ld 3/2 ) 



neutrons (ld 5/ 2-ld3/2) 



SkP 



protons (ld 5/ 2- 10*3/2) 



III li III Ml j ill a il 



CjoCfoCjoCfo C C* Co "Cf "C™ C„ c\ C" m C' w C' m C", B C'„ C\ Cj'Cf'Cj'cy 1 Cj Ci 

Coupling constants 

FIG. 13: (Color online) Same as in Fig. [6] but for the SkP 
Skyrme functional. 
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FIG. 14: (Color online) Same as in Fig. [g] but for the SkO' 
Skyrme functional. 



functionals, respectively. Comparison with the SLy5 re- 
sults previously presented in Fig. [HI shows that for all the 
three forces, patterns of regressions coefficients are very 
similar indeed. Small differences in the overall magni- 
tudes are probably related to differences in the power a 
of the density dependence, Eq. ((U, which is equal to 1/6 
(1/4) for SLy5 and SkP (SkO') and in the values of the 
effective mass, m*/m = 0.7, 0.9, and 1 for SLy5, SkO', 
and SkP, respectively. 



C. Fits to experimental data 

The regression coefficients Ii m obtained in the present 
work can be used to evaluate the possibility of readjusting 
the EDF coupling constants so as to better describe the 
experimental data. Indeed, if for a given set of coupling 
constants C m the self-consistent s.p. energies are equal to 
Ci, then corrections to coupling constants AC m = C** — 
C m are given by: 



cxp 



(12) 



The matrix 2j m being rectangular, it cannot be simply 
inverted, but, instead, one can use its singular value de- 
composition [1 81 ] to determine the result, i.e., 



w n V„ 



(13) 



where columns of matrices U and V are orthogonal and 
singular values w n are positive. In our case, the num- 
ber of experimental data q xp is much larger than the 
number of coupling constants C m , which is equal to 12. 
Therefore, the possibility of describing experimental data 
clearly depends on whether the column Aej can be ex- 
pressed as a linear combination of 12 orthogonal columns 
of matrix U. In any case, the result of the least-square 
fit to experimental data is given by 



AC' 



(14) 



for n max =12. 

Singular values w n obtained for regression coefficients 
Iim corresponding to three Skyrme functionals arc shown 
in the lower panel of Fig. 1151 Independently of which 
functional is concerned, the singular values decrease ex- 
ponentially from about 1.3 to 0.001. Again, this illus- 
trates a great universality of the regression coefficients. 

On the one hand, as seen in Eq. ([13]) , small singular 
values w n correspond to columns of matrices U and V 
that have very small influence on the matrix of regression 
coefficients Ii m . On the other hand, these same small 
singular values have very large impact on the coupling 
constants in Eq. (fT4")l . Therefore, expression (fF4"|) should 
be used for non-maximum values of n max , i.e., by cutting 
off the small singular values. By that, one removes from 
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0,001 

2 4 6 8 10 
Number of singular value 

FIG. 15: (Color online) Bottom: singular values correspond- 
ing to the regression coefficients h m calculated for the SkO', 
SkP, and SLy5 Skyrme lunctionals. Top: the rms deviations 
of the s.p. energies from experimental data, obtained for a 
given number of singular values kept in Eq. (|14[) . 



the final result those features of the regression coefficients 
that are poorly defined, and at the same time one keeps 
only those linear combinations of coupling constants that 
are well determined by data. 

In our analysis, for experimental data we adopted 
M = 58 proton and neutron negative s.p. listed in 
Ref. [l9j]. Since here we study bare theoretical s.p. en- 
ergies, we have also subtracted from the experimental 
values the calculated mass polarization corrections re- 
lated to the one-body center-of-mass correction [5j. In 
the upper panel of Fig. 1151 we show the rms deviations, 




1/2 



(15) 



of the fitted s.p. energies, obtained for n max = 
1,2, ... , 12. For the standard parametrizations studied 
here, one obtains Ae rms = 1.73, 1.61, and 1.45 MeV 
for SLy5, SkO', and SkP functionals, respectively. At 
?^max = 1, one obtains the rms deviations that are only 
marginally better. With increasing n max , the results 
gradually improve, and at n max = 12 (for the full least- 
square fit) they reach limiting values of about 1.1 MeV. 

Although in particular cases, some combinations of 
the coupling constants are especially important, i.e., at 
?^max = 7-8 for SkO' or at n max = 4 for SkP, the rms 
deviations decrease rather steadily, and do not show pro- 
nounced affects at low values of n max , observed previ- 
ously for the total binding energies, see Fig. 2 in Ref. [9j . 
On the one hand, this means that the experimental val- 
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FIG. 16: (Color online) Histograms of residuals, Eq. (|12|l . for 
standard (upper panels) and refitted (lower panels) Skyrme 
functionals SkO', SkP, and SLy5. 



ues of the s.p. energies may, in principle, determine val- 
ues of more coupling constants than those of the total 
binding energies. On the other hand, the limiting value 
of Ae rms = 1.1 MeV, which can be obtained by refitting 
the standard functionals, is rather disappointing. Fig. 1161 
shows histograms of residuals, Eq. (fl"2"|) . for standard (up- 
per panels) and refitted (lower panels) Skyrme function- 
als. One can see again that the improvement obtained 
by the fit is not very spectacular, with quite a number 
of s.p. energies differing from experiment by 1 MeV or 
more. This indicates, that to obtain a better agreement 
with experimental one may have to extend the form of 
the standard Skyrme functional. 



IV. CONCLUSIONS 

The present work constitutes a step towards the search 
for a spectroscopic-quality energy density functional, 
which would properly and precisely account for the 
single-particle structure of atomic nuclei. We analyzed 
dependence of single-particle energies on coupling con- 
stants of the Skyrme functional and showed that this de- 
pendence is almost linear. Such an observation allowed 
us to focus on linear regression coefficients, which charac- 
terize variations of individual single-particle states with 
changing coupling constants. 

We discussed the regression coefficients obtained for 
single-particle states in 9 doubly magic spherical nuclei. 
The results turn out to be rather generic, with only a 
weak dependence on the particular variant of the Skyrme- 
functional parametrization. We showed that for nuclei 
near stability, isovector coupling constants have signif- 
icantly smaller impact on single-particle energies than 
isoscalar ones. We also showed that the effective-mass 
coupling constant Cq cannot be considered as merely 
changing the overall density of states, which was up to 
now a commonly accepted view. This coupling constants 
also significantly influences relative positions of single- 
particle levels, including the splitting of spin-orbit part- 
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ners. We explained these effects by a detailed analysis 
of densities and mean fields related to this coupling con- 
stant. 

By applying the singular value decomposition to the 
matrix of the regression coefficients, we performed fits of 
calculated single-particle energies to experimental data. 
We showed that even by fitting all the 12 coupling con- 
stants of the standard Skyrme functional to experiment, 



one is unable to obtain the rms deviations of single- 
particle energies below about 1.1 MeV. This discouraging 
result points out to a necessity of extending the form of 
the Skyrme functional beyond the standard parametriza- 
tion. 

This work was supported in part by the Academy of 
Finland and University of Jyvaskyla within the FIDIPRO 
programme and by the Polish Ministry of Science. 
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